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Abstract. Let z\,...,Zk be distinct grid points. If fk,o is the prescribed value of a 
function at the grid point Zk, and fk, r the prescribed value of the r-th derivative, for 
1 < r < n>k — 1, the Hcrmitc intcrpolant is the unique polynomial of degree N — 1 
(N = m + • • • + Tlx) which interpolates the prescribed function values and function 
derivatives. We obtain another derivation of a method for Hermite interpolation re- 
cently proposed by Butcher et al. [Numerical Algorithms, vol. 56 (2011), p. 319-347]. 
One advantage of our derivation is that it leads to an efficient method for updating 
the barycentric weights. If an additional derivative is prescribed at one of the inter- 
polation points, we show how to update the barycentric coefficients using only O (N) 
operations. Even in the context of confluent Newton series, a comparably efficient and 
general method to update the coefficients appears not to be known. If the method is 
properly implemented, it computes the barycentric weights with fewer operations than 
other methods and has very good numerical stability even when derivatives of high 
order are involved. We give a partial explanation of its numerical stability. 



I. Introduction 

An example of a Hermite interpolant in barycentric form is the unique polynomial 
ir(z) of degree 3 such that 7r(— 1) = 7r'(— 1) = f_ ± , 7r(l) = f\, and 7r'(l) = f 1 as 
given by 
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In general Hermite interpolation, the function value and its first — 1 derivatives are 
prescribed as 

fk,0i ■ ■ ■ i fk,n k -l 

at the interpolation point or grid point Zk for each Zk from the list z±,...,Zk- The 
problem is to find a polynomial ir(z) of degree N — 1, where N — n\ + • • • + uk, such 
that 



(1.2) 



d r 7i(z) 



dz r 



fk, r for r = 0, . . . , n k - 1 



z=z k 
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at each grid point z&. We shall always assume the grid points Zk to be distinct. For an 
elegant proof of the existence and uniqueness of the Hermite interpolant, see [3]. 

The polynomial 7r(z) can be represented in either the Newton form or the barycentric 
form. In the Newton form, the grid points Zk must be ordered in some way [2] . If the 
grid points are not carefully ordered, the Newton form is susceptible to catastrophic 
numerical instability (HE]. In contrast, the barycentric form does not require the grid 
points to be ordered and treats all the grid points equally. 

Barycentric form. The barycentric form of the interpolant is 



tt(z) = n*(z) f; fk > n *-\ ( Wk >° ) + Z^- 2 ( Wk >° ^ + W ^ 



1.3) 



T[ (n k - 1)! \(z - z k )) (n k - 2)! \(z - z k ) 2 (z 

w k ,o 




Zk) nk 




where tt*(z) is defined as Ylk=i{ z — Zk) nk . It is evident from inspection that 7r(z), as 



represented in (1.3), is a polynomial of degree N — 1. For a unique choice of weights 



Wk, r , 7r(z) will satisfy the interpolation conditions (1.2). Determining and updating the 



weights (or coefficients) Wk, r of the barycentric form is the topic of this paper. 



The representation of 7r(z) shown in (1.3) is usually termed the first barycentric form. 



For the second barycentric form, take fk$ = 1 and f^ r = for r > 1 corresponding to the 
function f(z) = 1. By the existence and uniquen ess o f the Hermite interpolant, we have 



1 = tt*(z) J2k=i Sr^o 1 w k,r(z — z k ) nk r . Dividing (1.3) by the barycentric representation 
of 1, we have 

(-t a\ ( \ _ ^k=l 2^s=0 gl Z^r=0 UJk,r\Z Zk) 

This is the second barycentric form. The second barycentric form has a useful property 



that we now state. If the first barycentric form (1.3) satisfies the interpolation conditions 



(1.2), so does the second barycentric form (1.4). Less obviously, for any choice of the 



weights Wk, r with Wk,o ^ for k — 1, . . . , K, the second barycentric form is a rational 



function of z which satisfies the interpolation conditions (1.2) (see [S]). The second 
barycentric form is more robust in the presence of rounding errors in the weights Wk, r > 
as one may expect and as we will demonstrate in Section 4. 



For deriving the weights Wk, r , we shall work with the first barycentric form (1.3). 
It is obvious from inspection that either of the two forms can be used to evaluate the 
interpolant n(z) at a given point z using O(N) arithmetic operations, although the 
second barycentric form is more robust in the presence of rounding errors. 

Calculating barycentric weights using divided differences. One of the methods 
for finding the weights Wk )T in 0(N 2 ) operations is due to Schneider and Werner [8]. Two 
ideas go into that method. We briefly describe the two ideas in the simpler context of 
Lagrange interpolation (see Werner [T3] for this special case) where the problem is to 
find a polynomial p(z) of degree K — 1 such that p(z k ) = fk,o holds at each grid point 
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The Newton representation of p(z) is ao + a%(z — Z\) + a<i{z — Z\){z — z<i) + • • • + 
a K-i I\k<K( z ~ z k) for some coefficients a^. The Lagrange representation takes the form 



z — z 



3) 



J2k=i bkLk(z), where Lk{z) is the un- normalized Lagrange cardinal function Tlj^ 
The first idea is to note that the coefficients and the coefficients bk are connected by a 
triangular matrix [5] . The second idea is to note that if p(z) = 1 then ao = 1 and = 
for k > 1 and the coefficients bk are equal to the Lagrange weights Wk- The triangular 
relationship between the two sets of coefficients is inverted using divided differences to 
compute the Lagrange weights. The method of Schneider and Werner for computing the 
barycentric weights Wk, r of the Hermite interpolant follows the same logic. 

We make two comments about the method of Schneider and Werner. Firstly, because 
the Newton expansion is used implicitly, the method depends upon divided differences 
and finding a good ordering of the grid points z^. Even if the grid points are carefully 
ordered and scaled using the logarithmic capacity of the interval of interpolation, it is 
unclear from the extant literature if the method is numerically stable in the presence 
of high order derivatives. As our discussion in Sections 4 and 5 will indicate, numerical 
stability is a more delicate issue when high order derivatives occur in the Hermite data. 

Secondly, the method for finding barycentric weights described in Section 2 typically 
has a lower operation count than the method of divided differences. However, Schneider 
and Werner [S] allow a general denominator in the barycentric form while we specialize 
the denominator to be 1. Even though it could be possible to specialize the method 
of divided differences to lower the operation count, that is unlikely to make it a better 
method for computing the barycentric weights. The argument for the method of Section 
2 is its simplicity, directness, and numerical stability. 

The method of Butcher et al. [1J and its extension. In Section 2, we give another 
derivation of the method of Butcher et al. for computing the barycentric weights. The 
method was derived by Butcher et al. using contour integrals and the manipulation of 
infinite series. Our derivation is more direct. The formalism of Butcher et al. extends 
to Birkhoff interpolation (or Hermite interpolation with incomplete data), a problem 
which is not discussed here. 

To explain the basic idea we use in Section 2 for computing the barycentric weights, 
we go back to the Hermite interpolant shown in (1.1). We want to compute a cubic 
polynomial whose Taylor expansion to two terms with center z — 1 is equal to fi+fx(z — 
1) and with center z = — 1 is equal to f_i + f_ x {z + 1). In the prefactor (z — l) 2 (z + l) 2 , 
(z — l) 2 has an effect on the Taylor expansion at z — 1 that is easily neutralized by 
dividing by (z — l) 2 . The Taylor expansion of (z + 1)~ 2 to two terms about z — 1 is 
1/4 — (z — l)/4. Therefore if (z + l) 2 is multiplied with 1/4 — (z — l)/4, the product is 
1 + O ((z - l) 2 ) . Thus the product 



(* + l ) 2 ( z ~ >< (i~Ti)2 X (i " V) X (/l + f[{z ' l)) 

- <* + ^ - ^ (* (4(^17 - 4(^1) ) + 4A) ) + «* - ^) 
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is equal to fx + f[(z — 1) + O {{z — l) 2 ). This observation along with the fact that the 
expression above is evidently O {{z + l) 2 ) explains part of the barycentric representation 
shown in (1.1). 

In general, computing barycentric weights comes down to finding the coefficients in 
the Taylor polynomial of expressions of the form Y[^ k {z + z\. — Zj)~ nj . The coefficients 
of such expressions can be found efficiently and with good numerical stability using a 
standard identity from symmetric function theory. 

In Section 3, we show how to update the barycentric weights when an additional 
data item is introduced using only O (N) operations. Such an efficient update is not 
available for the Newton representation of the Hermite interpolant [2J. If the Newton 
representation orders the grid points as z 1: ■ ■ ■ , Zk, it can be easily updated if a function 
value is introduced at a new grid point z K+1 or if an additional derivative is introduced at 
Zr- However, introducing an additional derivative at one of the points Zx, ■ ■ ■ , zk-\ will 
require a lot of the divided differences table to be recomputed. Indeed, if the additional 
derivative is introduced at zx, the entire divided differences table has to be recomputed. 
The method in Section 3 handles all these cases using only O(N) operations. 

Numerical stability. The second barycentric interpolant (1.4) of the Runge function 
l/(l+z 2 ) in the interval — 1 < z < 1 using K = 512 grid points and n—1 = 47 derivatives 
has a numerical error of about 10~ 15 in the oo-norm and in double precision arithmetic. 
If the first barycentric form is used, the error is less than 10 -12 everywhere except very 
close to the endpoints. The first barycentric from is very sensitive near the endpoints 
z = ±1. At the endpoints, the error in the first barycentric form is as large as 10 -7 , a 
phenomenon that is discussed in Section 4. Here we mention that the endpoints ±1 are 
not included among the 512 grid points. 

It is probably unwise to restrict the discussion of numerical stability to the interpo- 
lation error ||7r(z) — /(z)||oo- This is particularly true for the second barycentric form 
which is quite good at producing an accurate interpolant even if the barycentric weights 
Wk, r themselves are inaccurate. In Section 4, we use extended precision computations to 
check the relative errors in the barycentric weights. We find that the barycentric weights 
are computed with good accuracy. 

The key to the accuracy of the barycentric weights appears to be the behavior of a 
triangular system that appears in the method of Section 2 (see Lemma [2] in particular). 
In Section 5, we prove that the barycentric weights are computed with small relative 
error in a few limited situations. 



2. Barycentric weights 
The Hermite interpolant tc(z) is the unique polynomial of degree N—1 which satisfies 



the N — 1 interpolation conditions (1.2). The number of interpolation conditions at the 
grid point z^, 1 < k < K, is and N = rix + ■ ■ ■ + rajc- The polynomial n*(z) defined 
as Y\k = x{z — Zk) nk is of degree N. The Hermite interpolation conditions can be reworded 
as requiring the Taylor expansion of ir(z) centered at z = z k to be equal to 

(2.1) fk,o + / M (* "**) + •••+ (^Ti)! ^ ~ Zk ^ +°(( z - ^) nfe ) 
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at each grid point z k . 

Define ir k (z) = tt*(z)(z — z k )~ nk so that 7t k (z) is a polynomial of degree iV — n^.The 
polynomial 



r=ni— 1 



(2.2) W k (z)= w k , r (z-z k y 



r=0 



is defined by the requirement ir k (z)W k (z) = 1 + O ((z — z k ) nk ) . In other words, W k (z) 
is equal to the Taylor series of 7[ k (z)~ l with center z = z k and truncated at the n^-th 
power. 

The polynomial 

rr k (z)W k (z) U + f kA (z - z k ) + ■ ■ ■ + £^y ( z ~ **) n * _1 ) 



has a Taylor expansion at z — z k which is identical to (2.1) because 7t k (z)W k (z) — 1 + 
O ((z — z k ) Hk ) in the limit z — > z k . In addition, for j ^ k the polynomial is O ((z — Zj) nj ) 
in the limit z — > Zj because n k (z) is O ((z — Zj) n: >) in that limit. However, the degree 
of the polynomial is N + n k — 2 which is more than N — 1 if n k > 1. This difficulty is 
easily fixed. We simply need to multiply W k (z), which is thought of as a polynomial in 
[z — z k ), and Y^r^o 1 ^~( z ~ z k) r and drop all terms of order (z — z k ) Uk or higher. The 
following polynomial 

(2.3) 7t k ( Z ) (f kfi (w kt0 + ■■■+ W Knk -i{z - Z k ) nk ^) +■■■ 

. fk,n k -2(z - Zj,)"*^ 1 fk,n k -l(z ~ Z k ) 7lk ^ 1 

+ ( w k,0 + w^z - z k )) + w kfi 

is of degree N — 1. For j ^ k, it is O ((z — Zj) nj ) near Zj because ir k (z) is divisible by 
(z — Zj) nj . At z — z k , it satisfies the interpolation condition (2.1). We have the following 
lemma. 

Lemma 1. The barycentric weights w kr , r = 0, 1, . . . , n k — 1, are the coefficients of the 



unique polynomial W k (z) of degree n k — 1, shown in (2.2), which satisfies ir k (z)W k (z) 
1 + O ((z — z k ) nk ) in the limit z — )■ z k . 

To calculate the barycentric weights, it is useful to look at the following equation 

(2.4) 7 V7 ^ v = l+X lZ + X 2 z 2 + 

(1-^)...(1-^ X 

aij V "2/ V a N 



We assume a k ^ andXi,X 2 , ... on the right hand side of (2.4) are defined by expanding 
the left hand side in powers of z. 
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To find an efficient method to compute the X r , we define V r = J2k=i a k T an d differen- 



tiate (2.4) with respect to z to get 



1 N 1 / z 

]T— 1 = X 1 + 2X 2 z + 3X 3 z 2 + ■ 



(l + X^ + X 2 z 2 + ■■■){V 1 +V 2 z + ■■■) = Xi + 2X 2 ^ + ■ 
Equating coefficients of 2, we have the following lemma. 



Lemma 2. The quantities X r defined by (2.4) are related to the inverse power sums 
V r = Eti < by 

Tx = V l 

2X 2 = V 2 +X{Pi 

3X 3 = V 3 +X l V 2 +X 2 V 1 

4X 4 = V^+X{P 3 +X 2 V 2 +X 3 Vi 

and so on. 



Coincidentally, Lemma [2] appears as Lemma 2 in pQ and in [TU] , but with different 
proofs. It has been well known in symmetric function theory for a very long time. 
If we note that 

z — z k \ 




*k{z) 1 = nizk-zj) _ 

Lemmas [T] and |2] imply an algorithm for computing Wk, r - The algorithm is to begin by 
computing Ck = I\j^k( z k — z j)~ nj and define V r = J2j^k n j( z j — Zk)~ r ■ The quantities 
X r are obtained from Lemma [2] and Wk, r = CkX r . The grid points Zk are assumed to be 
distinct. 

Pseudo-code derived from a C++ implementation is displayed as Algorithm [TJ For 
another derivation of the same algorithm, see Butcher et al. pQ. 

The functions PowerSums and InversePoly defined on lines 7 and 15 of Algo- 
rithm [l] perform 2nK and n 2 arithmetic operations, respectively. Summing the cost of 
invoking those functions from lines 30 and 31, we find that the number of arithmetic 
operations performed by Algorithm [l] is 2NK + J2k=i n k to leading order. Roughly half 
of these operations are additions or subtractions and roughly half are multiplications. 
On line 29, a total of K 2 divisions are performed. That count can be easily reduced 
to K 2 /2 because the differences Zk — Zj that needed to be inverted are only K 2 /2 in 
number. The arithmetic operations on lines 32 and 33 do not affect the leading order of 
the total count. 

The method of divided differences for computing the barycentric weights Wk, r requires 

K{K — l)/2 divisions and ^ k multiplications [8]. The number of subtractions or 

additions is about the same as the number of multiplications. To compare the operation 
counts of the two methods, consider the situation in which n\ = ■ ■ ■ = = n. In that 
situation, the number of multiplications used by the method of divided differences is 
(n 2 K 2 — Kn 2 ) /2 as against nK 2 + n 2 K used by the method of Newton identities. The 
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Algorithm 1 Computing barycentric weights 

1: function findPI(£i, . . . , Ck, ni, . . . , Uk) 
2: return d" 1 • • • Ck 

3: end function 

4: function nInvertLiSt(^i, . . . , z K , (i, • • • > Ck) 
5: ( k = -l/z k for k = 1, . . . ,K 

6: end function 

7: function PowerSums(0, . . . , Ck, n-i, • • • , n K , Pi, ... , V n -\) 
8: Temporaries: ti,...,tx 
9: if. = n k for 1 < k < K 
10: for r = 1, . . . , n — 1 do 

11: 4 = ( fc t fc for 1 < k < K 

12: V r = h + ■ ■ ■ + t K 

13: end for 

14: end function 

15: function InversePoly(P 1; . . . ,V n -i, X ,Xi, . . . ,X„_i) 

16: X = 1 

17: for r = 1, . . . , n — 1 do 
18: X r = ELi %- s /r 

19: end for 

20: end function 

21: Comment: The grid points Zfc must be distinct. 

22: function HermiteWeights(zi, . . . , z K ,ni, . . . w k , T with 1 < k < K and < 

r < n k ) 
23: for k = 1, . . . , if do 

24: Temporaries: z£> ■ ■ ■ > z^-i* n'^ . . . , n'^^n, C k , Ci, • • • , Gf-i 

25: n = n k 

26: Temporaries: Pi, . . . , TV-i) X , . . . ,X n _i 

27: z'a = z k — Zj for 1 < j < k and Zj = z k — Zj + i for k < j < K — 1 

28: n'j = rij for 1 < j < k and n'j = rij + \ fork<j<K — l 

29: NlNVERTLlST(^, . . . , Z^-nCl, ■ ■ • , C*-l) 

30: POWERSUMS(Ci, • • • , Ctf-lX, • • • , "jr-i, • • • > ^n-l) 

31: InversePoly(Pi, . . . ,V n -l, X ,Xi, . . . ,X n _i) 

32: C k = 1/findPI(^, . . . , • • • > "k-l) 

33: U7fc,r = C k T r for r = 0, 1, . . . , n — 1 

34: end for 

35: end function 



ratio of the two quantities is (K — 1)/ (2K/n+ 1). The same ratio applies to the number 
of subtractions or additions. The two methods have similar operation counts if n = 2, 
but for n > 2 Algorithm [T] uses fewer arithmetic operations by a factor of approximately 
n/2. 
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The argument for Algorithm [T] over the method of divided differences is its simplicity 
and numerical stability, which will be illustrated in Sections 4 and 5. In addition, 
Algorithm [T] has a favorable operation count. 

3. Updating the barycentric weights 

Suppose the barycentric weights w k>r have been computed assuming that the function 
value and n k — 1 derivatives are prescribed at grid points z k , for 1 < k < K. Suppose a 
new item of data is introduced at the grid point (. If ( is a new grid point, the new item 
of data is the function value at the new grid point. If £ = z k for some 1 < k < K, the 
new item of data must be the n k -th derivative at z k . In either case, one new barycentric 
weight must be computed and all the other barycentric weights must be updated. We 
show how to do that in O(N) operations, where N = ni + ■ ■ ■ + nx- 

Suppose new data is introduced at ( and Zk ^ (■ Updating the weights Wk, r , < 
r < n k , is the easy part. The polynomial Tr k (z) must be changed to n k (z)(z — (). The 
weights Wk, r were determined such that n k (z)W k (z) = 1 + O ((z — Zk) nk ) for z — > z k 
with Wk(z) = Z)r=o w k,r{z — Zk) T '. If the new weights are w' kr and the corresponding 
polynomial is W k (z), then we require 

n k (z)(z - C)W h (z) = 1 + O {{z - ZkD . 
To ensure this condition, it is enough if we find W' k (z) of degree n k — 1 such that 

(z - C)W' k {z) = ((z - z k ) + (z k - Q)W k {z) = W k (z) + 0((z- Zk) n *) 
in the limit z — > Zk- This gives the equations 

(z k - ()w' kfi = w kfl 
(3.1) (z k - ()w' kr + Wfc ir _! = w k , r for 1 < r < n k 

which are solved to deduce the updated barycentric weights w' k r . It takes two operations 
to update each barycentric weight w kyr for k with Zk ^ (■ 

If z k ( for 1 < k < K, then ( is the new grid point z^+i- The weights w ktT are 



updated using (3.1 ) for 1 < k < K as already described. We need to compute the new 



weight wk+i,o- This weight is computed using the formula wk+i,o = Uj=i(C ~ z k)~ Uk - 
If £ is a new grid point, namely zk+i, the total cost for updating all the weights and 
computing the new weight is O(N). 

If ( = z K for some k, 1 < k < K, the weights w kyT are updated as already described 
if k 7^ k. We are left with the case of updating the weights w K>T with ( = z K for some 
K, 1 < k < K. In this case, the new data item is another derivative at z K . Efficient 
updating of the weights w KtT causes complications that we now turn to. 

Define 

A k = ll/(zj — Zk) repeated rij times j ^ k and 1 < j < 
C k = \{{z k - Zj )- n K 

Here Ak is a multiset. Define V r (Ak) as the sum of the r-th powers of the elements of Ak- 
Let X r (A k ) be related to V r (A k ) by the triangular identities of Lemma 12} According to 
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Algorithm [TJ the weight w k)r is given by w k>r = C k X r (A k ) for < r < n k and 1 < h < K. 
When computing the weights w kj0 , ■ ■ ■ ,Wk, nk -i, the 2n k — 1 intermediate quantities 

(3.2) C k ,Vt(A k ), . . .,V nk -x(A k ),li(A. k ), . . . ,l nk -i(A k ) 

arise. We assume that these intermediate quantities are stored for each k. The total 
number of intermediate quantities stored is 2N — K. 

The reason for storing these intermediate quantities is as follows. Suppose a new 
data entry is introduced at the grid point £ = z K . Then a new weight w K) n K needs to 
be generated. To calculate that weight, I nK (A K ) will be generated using the triangular 
identities of Lemma [2] (the update rules of ( 3.1| are of no use here). The intermediate 



quantities are used to calculate I nK (Ac)- 

When new data is introduced at the point (, the intermediate quantities are updated 
at all grid points z k regardless of whether z k ^ £ or z k — C (so than we can handle a 
new data item introduced at z k later on). If z k ^ (, some of the intermediate quantities 



in the list (3.2) are updated as follows. 

C' k <- C k /(z k -() 
(3.3) V r (A' k ) 4r- V r (A k ) + (C - z k y r for l<r<n k 

where the primes denote updated quantities with A' k = A k U {1/(C ~ z k)} (this is a 



multiset union not a set union). To update X r , we go back to ( |2.4 ), which defines 
Xi,X 2 , . . ., and deduce the following: 

(3.4) l r (A' k )- Ir - l{A ' k) =X r .(A k ) for r = 1, . . . ,n k - 1 

S — z k 



where it is assumed that X = 1. Using (3.3) and (3.4), it costs two operations to update 



each intermediate quantity in the list (3.2). Since the total number of intermediate 



quantities is O(N), the cost for updating the intermediate quantities for k ^ k is also 
0(N). 



In £ = z K , A K does not change and none of the intermediate quantities in (3.2) with 
k = K needs to be updated. However, we need to generate V nK (A K ) and I nK (A K ). These 
are generated using the formulas 

n K 

i¥=K 8=1 

the first of which is by definition and the second is a triangular identity of Lemma [2] If 
the temporaries tj that arise on line 8 of Algorithm [T] are preserved, V nK can be computed 
using only 2nn — 1 operations. The barycentric weights u> K ,o, • • • , w K ,n K -i do not change. 
The weight w Ki „ K is computed as C K I nK . We have thus completed describing an 0(N) 
method for updating all the barycentric weights and generating a new weight when a 
new data item is introduced. 

If the barycentric weights are generated by repeatedly adding new items of data at 
grid points, the grid points must be ordered in some fashion. As we noted in the intro- 
ductory section, ordering the grid points may make the method vulnerable to numerical 
instability. If the weights need to be updated a small number of times to accommodate 
new items of data, the method described in this section will be numerically stable. 
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A more general updating problem is to update the existing barycentric weights and 
generate new ones, when n' k new data items are introduced at the grid point z k for 
1 < k < K. Algorithm [T] and the updating algorithm of this section can be combined 
to deduce an efficient and numerically stable solution to this more general updating 
problem. 



4. Illustration of numerical stability 

If the grid points z k , 1 < k < K, and the number of derivatives at grid points n^, 
1 < k < K, are chosen without care, the barycentric interpolation of f(z) using the 
Hermite interpolant tt(z) will be very badly conditioned. For such grids, it makes little 
sense to speak of the numerical stability of Hermite interpolation using either the first or 



second barycentric forms given by (1.3) and (1.4), respectively. We must first ascertain 
that the choice of Zk, nk leads to a well-conditioned interpolation problem. 

Let lk, r (z) denote the polynomial of degree N — 1 (N = n\ + ■ ■ • + tik) such that 

d s lk,r{z) 



dz & 



= 8 riS 8k t i for 1 < / < K, < s < ni. 

Z = Zl 

The fundamental polynomials (or cardinal functions) of Hermite interpolation lk, r (z) 
have a prescribed behavior at the grid points. However, they may become very large 
in-between the grid points leading to a poorly conditioned interpolation problem. 

If the grid points are chosen as Zk = cos((2/c — V)7r/2K), 1 < k < K, (these are the 
Chebyshev points), and n\ = • • • = n K = n, Szabados [TU] has proved upper bounds on 
the fundamental polynomials in the interval [—1,1]. In particular, 

( (\ogK\ . 

-Si? Mz)| - r i \ 

h=1 [ O ij^:) if n — r is even 

for r = 0, . . . , n — 1. These bounds imply that a small change in the interpolation data 



fk, r will result in a small change in the interpolant 7r(z) (defined by (1.3) or (1.4)) thus 
showing the Hermite interpolation problem to be well-conditioned. 

In all our computations, the Chebyshev points z\. are replaced by 2z k and correspond- 
ingly the data fk, r are replaced by fk, r /2 r - If 2 £ [—1, 1], we have 

&(n <*-*)) = 5 

because the capacity of an interval is a quarter of its length [31 p. 83]. So we may expect 



the prefactor tt*(z) in the first barycentric form (1.3) as well as the quantity Ck, which 
occurs on line 32 of Algorithm [TJ to be roughly of the order l/2 nK . For reasonably large 
K and n, Ck and ir*(z) underflow in IEEE double precision arithmetic. Replacing z k 
by 2zk and taking z G [—2, 2] changes the limit above to 1 implying that the quantities 
such as C k that occur in the algorithm are better scaled. 

Barycentric Hermite interpolation problem is highly susceptible to overflows and un- 
derflows. Scaling the Chebyshev points so that the capacity of the interval is 1 as in the 
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Figure 4.1. (A): Interpolants tt(z) to the Runge function f(z) = 1/(1 + 
z 2 ) with K = 4 and n = 1,2,3. (B): Dependence of interpolation error 
on the number of grid points K for fixed number of derivatives n at each 
grid point. 



previous paragraph is only a partial cure. Even if the Zk are scaled as indicated, a prod- 
uct evaluated in the order (2z — 2zi), (2z — 2z\)(2z — 2z 2 ), (2z — 2z\)(2z — 2z 2 ){2z — 2z 3 ), 
and so on may overflow or underflow in its intermediate stages even though the final 
result can be represented in machine arithmetic. To decrease the chances of such a thing 
happening, we reorder the Chebyshev points Zk in the Leja ordering [7] in addition to 
multiplying them by 2. 

In our implementation of the barycentric formulas (1.3) and (1.4), the function values 
input to the interpolation procedures are fn,r/r\ not f n ^ r . This eliminates the need to 
first multiply and then divide a quantity by r! (the computation of which is obviously 
prone to overflows) for certain functions such as the Runge function discussed below. 

Typical C and C++ compilers implicitly turn on the flush-to- zero optimization. Given 
the tendency of Hermite interpolation to form intermediate quantities that are prone to 
underflow, this optimization must be explicitly turned off. Another compiler option that 
is implicitly used on current x86 platforms is to replace the usual IEEE division by a less 
precise but faster floating point division. Intel/ AMD manuals give no information about 
the loss of precision or the gain in speed. We have not encountered a single problem 
where this so-called "optimization" makes the computing faster (nor have we noticed a 
loss in precision). There is no point in using it. The Intel C/C++ compiler options that 
must be used are no-ftz and either prec-div or fp-model strict. 

We now turn to the Runge function f(z) = 1/(1 + z 2 ), which is the subject of com- 



putations shown in Figure 4.1 A numerically stable method to calculate higher order 



BARYCENTRIC HERMITE INTERPOLATION 



12 




-1 -0.5 0.5 1 10° io' io z io 3 io* 



z N=nK 

(A) (B) 

FIGURE 4.2. A): Interpolants n(z) to the hat function f(z) = 1 — \z\ 
with K = 4 and n — 1, 2, 3. (B): Dependence of interpolation error on the 
number of grid points K for fixed number of derivatives n at each grid 
point. 



derivatives of the Runge function for z G [—1,1] is implied by the following calculation: 

/(*) = 
d r f(z) 




dz r 2i\(z-i) r+l (z + i) r+1 ) 

If z — % — Rexp(i9), then we have f( r '(z)/r\ = (— l) r ~ 1 i? r+1 sin(r + 1)6. This formula 
has excellent numerical stability for z G [—1,1]. 



Figure 4.1 shows that the error in interpolating the Runge function is c n exp(— anK). 
The positive constant a appears to be independent of n. However, c n increases with n. 
Therefore if we are allowed to use N items of information about the function to form an 
accurate interpolant, the best choice for this example is to take them all to be function 
values. 

For the hat function f(z) — 1 — \z\, Figure 4.2 shows that the error is proportional to 
1/K. If the error is represented in the form c n /nK, c n once again increases with n. 

The Runge function is analytic in an open set around [—1, 1], while the hat function is 
of bounded variation but not even continuously differentiable. In the case n = 1, which 
corresponds to ordinary polynomial interpolation, the differing rates of convergence are 
easily explained by the difference in smoothness and analyticity [3]. That connection 
likely persists for n > 1, which corresponds to Hermite interpolation. 



The interpolations discussed so far used the second barycentric form. Figure 4.3 



compares the first barycentric form (1.3) to the second (1.4). Even with K = 512 grid 
points and n = 48 derivatives there is no sign of numerical instability at all. The first 
barycentric form is less accurate but has small errors for z G (—1, 1). At z — ±1, the 
first barycentric form produces much larger relative errors of about 1CT 7 . It is well 
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Figure 4.3. Interpolation errors to the Runge function f(z) = 1/ (1 + z 2 ) 
with K = 512 grid points and n = 48 derivatives at each grid point. The 
first barycentric form is notably inaccurate at z — ±1. 



known that the Chebyshev polynomials are bounded by 1 for z G [—1, 1] but increase 
rapidly outside the interval. That phenomenon is magnified by several orders for the 
fundamental polynomials of Hermite interpolation. Thus we should expect the Hermite 
interpolant with n = 48 to be very sensitive near the endpoints of the interval. 

Another phenomenon noticeable from Figure 4.3| is that the errors are somewhat 
elevated near the center of the interval. The second barycentric form exhibits neither 
elevated error near the middle of the interval nor sensitivity at the end points. It has 
excellent numerical stability. 



In Figure 4.4, we directly examine the weights used for barycentric interpolation. 
These weights are computed using Algorithm [T] Figure 4A uses K = 16 and n = 16. 
The Chebyshev points Zk were multiplied by 2 but the display in the figure does not 
use Leja ordering. From part (A) of the figure, we see that the weights Wk, r are quite 
small especially when r = 0. Thus in spite of using an interval of capacity 1, the Ck that 
occur on line 32 of Algorithm [T] are getting quite small. Using an interval of capacity 
1 still allows considerable fluctuations in Ylk =1 (z — Zk) around 1 because it is only the 
1/K-tYi power which is guaranteed to converge to 1. In computing the prefactor Ck, 
these fluctuations are raised to the power n, with n being 16 in this instance. 



From part (B) of Figure 4.4, we see that the relative errors in the weights Wk, r are 
quite small. However, the relative errors are the highest for r = n and Zk near the middle 
of the interval. These phenomena will be partially explained in the next section. 

In both Figures 4A and 4J3, the relative errors in the barycentric weights were ob- 
tained by comparing double precision results of a C++ program with extended precision 
computations in Maple. Figure 4J3 shows the increase in the maximum relative error in 
the barycentric weights Wk, r with n, which is the number of derivatives used at each grid 
point. The increase is mild. Indeed, in part (B) of the figure, which uses K = 512 the 
increase is only linear. Algorithm [T] for computing barycentric weights appears to have 
very good numerical stability. 
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Figure 4.4. (A) Magnitude of the barycentric weights w^ r and (B) Rel- 
ative error in the weights (the maximum relative error is 2.86 x 10~ 12 . 
Both plots use K = 16 grid points and n = 16 derivatives at each grid 
point. 




5. Partial explanation of numerical stability 

On line 32 of Algorithm [TJ the weight w^ r is computed using w^ r = C^L r . The Ck 
are merely products of differences and will be computed with excellent relative accuracy. 
We will restrict ourselves to the error accumulation in the computation of X r . 
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We assume that {cti, . . . , a^} is a multiset of numbers and that P r = J2i=i a 7 r - The 
quantities I r are defined by (2.4) and computed using P r and the triangular identities 
of Lemma [2j For computing Wk ;r , the multiset is taken to be z\ — Zk repeated n x times, 
Z2 — Zk repeated n 2 times, and so on — see line 27 of Algorithm [Tj 



-53 



each 



If u is the unit round-off (u = 2 for IEEE double precision arithmetic), 
arithmetic operation of the type a + b will evaluate to (a + b)(l + 8) with \S\ < u in 
floating point arithmetic. If l + 9 n — T\i = i{l + 8i) pi , where \5i\ < u and pi = ±1, we have 
&n < 1m where 7„ is defined as nu/(l — nu) (this is Lemma 3.1 of [5] for whose validity 
it is assumed that nu < 1). This lemma is convenient for accumulating the effect of 
several rounding operations. 

Another convenience is the counter notation [SJ [9]. In this notation, (n) stands for a 
product of the form n?=i(l + ^«) P S where \Si\ < u and pi = ±1. Evidently, if (n) — 1 + 5 
then \8\ < 7„ (assuming nu < 1). In floating point arithmetic, (a + b) + c evaluates to 
a (2) +6(2) +c(l). 

If Pi = l/«i + -- - + l/ajv is computed in the order it is written, it evaluates to 
l/«i (AT) + l/a 2 (AO + l/a 3 (iV - 1) • • • + 1/«tv (2) . 

If a« is perturbed to «j (AT — z + 2), the computed result is the exact result. Thus the 
computation of V\ is backward stable. However, the computation of V\ can still have 



large relative errors if there are near cancellations. In Figure 4.5 (A), the relative error 
in V\ (the maximum over all z^ which occurs for z^ near the middle of the interval) 
is 1.18 x 1CT 11 . In part (B), which correspond to K = 512, the relative error in V\ 
is 6.42 x 10~ 9 . If there is a point at and the positive and negative grid points are 
symmetric, the exact value of V T is zero for r an odd number and the relative error is 
infinite. However, in this case, V r will have small relative error for even r and the only 
power sums that influence the X r correspond to even powers. If r is even, V r is always 
computed with excellent relative accuracy as there can be no cancellations. 

Since the error analysis of V r is relatively tractable and in view of the comments 
made in the previous paragraph, we will temporarily assume all the power sums V r 
to be computed with small relative errors. The more complex issue is the manner in 
which the errors are amplified when the quantities X r are computed using the triangular 
identities of Lemma [2] Figure 4J3 shows that the errors grow only very mildly when 
going from V r to I r . 

To understand the error growth, we recast the triangular identities of Lemma [2] into 
matrix form: 



(5.1) 



1 

-Pi 
-P 2 
-Vs 



1 

-Pi 
-P 2 



2 

-Pi 



V 



/ 



/ 1 \ 

V : J 



( 1 \ 






V : J 



If the condition number of the matrix in (|5.1|) is k, the norm- wise relative error in the 

2e/t 



X r is bounded by 



l- 



where e is of the order of u, the unit round-off. Unfortunately, 



this bound is useless because k can be very large. If a single oij is less than 1, which is 
always the case, the P r will increase exponentially with r. The condition number of the 
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system (5.1), retaining only n + 1 rows in that system, will increase exponentially with 
n. Perhaps componentwise bounds [5J, 7.2] can be more useful. 

Although k is quite bad, there is a simple way to get a matrix system equivalent to 



(5.1) with a small condition number. Suppose the a,s are all multiplied by a scaling 
factor s and changed to saj. Then V r and X T are transformed to P r = s~ r V r and 
I r = s~ r X r , respectively. If s is sufficiently large then \P r \ < 1 for 1 < r < n. 

The (n + 1) x (n + 1) matrix obtained by replacing the power sums V r by their scaled 
version P r has condition number bounded by n + 1. In general, triangular systems have 
condition numbers that increase exponentially even if all the off-diagonal entries are 
bounded by 1 in magnitude [12]. Here the situation is different because the diagonal 



entries increase in the order 1, 1, 2, 3, . . . We have the bound 

, 7 | = \P r + hPr-l + ■ ■ ■ + Ir-lPl\ < 1 + \h\ + ■ ■ ■ + \I r -l\ < j 



where the last inequality is by induction. Thus all entries in the inverse will be bounded 
byl. 

If a large scaling factor s is used to control the condition number of the matrix, the 
difficulty comes in when we use X n = s n I n to estimate the relative error in X n . The 
absolute errors in I r can be bounded by a quantity of the form n 5//2 w, or even nu if the 
scaling factor s is chosen to make the power sums small enough , where u is the unit 
round-off. However, if I n is much smaller than 1 in magnitude, the norm-wise bounds 
will imply a large relative error in I n and in X n . This difficulty can be removed if a scaling 
factor s can be found such that \P r \ < 1 for 1 < r < n and \I n \ is 0(1). The linear 



growth in relative errors in Figure 4.5 suggests that such a scale might exist. Proving 



the existence of such a scale will require estimates of the quantities V r and X T 



The discussion in this section partially backs up the finding of Figure |4.5| that the 
growth of relative errors in the barycentric weights with the order of the derivative is 
mild. The contribution to the errors made by V r is therefore quite significant. Even 
though the computation of each power sum V r is backward stable, the sum itself could 



be ill-conditioned. In fact, by looking up n = 1 in Figure 4.5, we see that the relative 
error in V\ is not so small. We implemented compensated summation in x86 assembly 
language and verified that it does not reduce forward errors. 

The following proposition allows us to verify the accuracy of some of the barycentric 
weights. 

Proposition 3. // the oti are all positive, the relative error in V r is bounded by 7 r +jv 
and the relative errors in X r are bounded by 73 r jv for sufficiently small unit round-off' 
and r < N. 

Proof. The computed quantity V r can be represented as 

a± r (N + r) + cq r (N + r - 1) + • ■ ■ + a^ r < r + 1 > . 



Because each «j is positive, the computed quantity can be represented as V r (N + r) 
thus proving one half of the proposition. 
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For the other half, we note that that X r = (V r + X{P r -i + • • • + X r ^{Pi) jr. By in- 
duction, the quantity computed for X r can be represented as 

V T (N + r) (r) + ■ ■ ■ + Tr-iPi (3(r - l)N) (N + 1) (r) 

r 

In this expression the (r) factors that accompany each term in the numerator are caused 
by the arithmetic operations for computing X r using X s for 1 < s < r — 1. The factors 
that accompany each term in the numerator can be multiplied and replaced by (3Nr). 
Since each term is positive, it follows that the computed quantity can be represented as 
X r (3Nr). □ 

This proposition implies that the weights Wk, r are computed with excellent accuracy 
when Zk is either the least or the greatest of the grid points. Indeed, we see from Figure 
4.4[ B) that the relative errors near the end points of the interval are very small. If 



Zk is the least of the grid points, each z'j , 1 < j < K — 1 that occurs on line 27 of 
Algorithm [I] is negative and non-zero. Since — 1/z'j corresponds to otj in Proposition 
[3j the proposition is immediately applicable. If Zk is the greatest of the grid points, 
each —1/z'j is negative and the proposition is not immediately applicable. However, it 
is evident from inspection that V r and X r are negative or positive according as r is odd 
or even. Therefore every sum in the triangular identities of Lemma [2] has terms that are 
all positive or all negative. Thus the bounds on relative errors proved in the proposition 
will apply. 
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